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ABSTRACT 

Middle atmospheric transient Maxwell current densities generated 
by lightning induced charge perturbations are investigated via a 

simulation of Maxwell’s equations. A time domain finite element 

analysis is employed for the simulations. The atmosphere is modeled as 
a region contained within a right circular cylinder with a height of 
110 km and radius of 80 km. A composite conductivity profile based on 
measured data is used with charge perturbations are centered about the 
vertical axis at altitudes of 6 and 10 km. The simulations indicate 

that the temporal structure of the Maxwell current density is 

relatively insensitive to altitude variation within the region 

considered. It is also shown that the electric field and Maxwell 

current density are not generally aligned. 
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INTRODUCTION 


The use of Maxwell currents and current densities to descrioe the 

electromagnetic response of the atmosphere is a fairly recent 

development. Krider and Musser and Nisbet have suggested that the 

thunderstorm is fundamentally a current source and should therefore be 

investigated in terms of current densities rather than electric fields. 

The current densities generated by charge perturbations associated with 

lightning consist of either conduction (resulting from charge movement) 

or displacement (cdE/dt) terms. The sum of these current densities are 

referred to as the Maxwell current density (J , 

J m). 

j = VxH = J + cd E/at (1) 

m p 


J = Maxwell current density 
m 

J = the sum of all conduction current 
P 

densities including source terms 
c<9E/3t = displacement current density 

The Maxwell current density has several properties that may be 
exploited to more accurately describe both the local and global effects 
of lightning on the atmosphere. Probably the most important of these is 
;hat the divergence of the Maxwell current density is zero (V # J m = 0) 
making it a solenoidai quantity. As a consequence of this, the lines of 
the Maxwell current density form closed loops. Along these Maxwell 
current density streamlines the displacement term is usually dominant at 
low altitudes with the conduction current density dominating the high 
altitude Maxwell current density. 

This solenoidai character of the Maxwell current density implies 
that the electrical parameters of the entire path of circulation 
strongly affect its response and conversely; if the Maxwell current 
density along the streamline is mathematically describable, the 
corresponding streamline’s electric field may be formulated by a simple 
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time domain integration. An analogy exemplifying this phenomenon can be 

made to a source driving a series ladder network with each element 

consisting of a parallel RC. The relative analogies in the circuit are 

as follows: J ~ I , voltage ~ E, CdV/dt ~ edE/dt, V/R ~ crE. It is 

m s 

obvious that the total impedance of the circuit limits the current 
response, analogous to the peak surge impedance along the streamline's 
path governing the behavior of the Maxwell current density. 
Therefore, if the limiting factors (peak surge impedance) could be 
identified and incorporated in an atmospheric electrodynamic model, the 
electric fields along the entire streamline could be easily derived. 

Modeling of the Thunderstorm 

The simulated thunderstorm’s electrical activity is sustained by a 

constant current generator that exists between the upper and lower 
charge centers (6 and 10 km respectively, upper center positive). 
Intra-cloud lightning, resulting from the accumulation of the generator 
charge, occurs at time intervals determined primarily by the charging 
current and amount of charge exchanged per flash. For the purposes of 
modeling, the effects of the constant current generator will be analyzed 
separately from that of the charge perturbation associated with 

intra-cloud lightning. The resulting steady state and transient 
solutions will then be superimposed to determine the total electrical 
response of interest [Baginski, 1987]. 

Several implicit assumptions used in the modeling should be noted. 

The net amount of charge accumulation at either of the charge centers due 
to the thundercloud generator current and the neutralization of charge 

by lightning was assumed to be zero. This required that the time 

average amount of charge neutralized by lightning equaled the amount 
deposited by the thundercloud current generator. In the modeling a one 
ampere charging current was used. Charge accumulation at the charge 
centers occurred until a breakdown field strength was obtained [pg 72, 
Uman, 1987]. The resulting steady state field mapping at the time 
breakdown was reached was used as the steady state component of the 
simulated field signatures. Specific details of the spatial and 

temporal distribution of the thundercloud's charge used in the modeling 
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The computer model used in the study was based on an adaptation of 
an earlier computer code that was used to simulate lightning induced 
field signatures in the middle atmosphere and ionosphere [Baginski et 
aL, 1988]. The major changes necessary were to the inclusion of the 

thundercloud’s current generator as an additional electrical source and 
use of a conductivity profile based on measurements. Since the earlier 
code simulated transient behavior that was observed at greater distances 
from the source and had much shorter durations, no changes were required 
in time step or spatial discretization used. 

Charging Mechanisms 

In thunderstorm research, the most difficult phenomena to explain 
have been the processes involved in cloud electrification. The 
difficulty is twofold: on the one hand, there are a large number of 

possible mechanisms responsible for charge separation and current 
generation; on the other hand, it is usually impossible to isolate such 
mechanisms and test each for its relative effect. Regardless of the 
mechanisms, what is known is that a thunderstorm is sustained by charge 
separation which can be approximated by net positive and negative charge 
centers. The height of the charge centers is somewhat affected by 
seasonal changes and the geographic location. Typical heights of i0 km 
for the upper charge center and 6 km for the lower charge center are 
wideiv found in the literature and were selected for this research [pp 
60-6 i, Oman, 19871. 

According to measurements, the amount of charge neutralized by 
ligntning, and the thunderstorm’s current generator will vary over a 
considerably larger range than the charge center’s heights. For 
instance, Kasimer u959j has measured values of cloud electrification 
currents from less than 0.1 ampere to 10 amperes. Lightning return 
stroke currents are reported to have an even larger range of values [pp 
161-162 Urn an, 1987b It was therefore necessary tc carefully consider 
the affects the temporal signature and magnitude of these sources would 
have on the simulations of interest prior to their specification in the 
model. 
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Since the research presented here focused on field recovery 
behavior at times ~ 100 ms after the cessation of the lightning current, 
the temporal variation of the lightning current was assumed to have 
little effect on the simulations of interest. In fact, after the 
modeling and comparative studies were completed, a simulation was done 
where the time variation of the charge neutralized (source of transient) 
by lightning was described mathematically by a simple impulse (S(t)) and 
resulted in field recovery curves that were identical. 

Sunde’s (1968) lightning current model was used in the study in the 
form shown below: 


I (t) * I0(exp(-at) - exp(-bt)) (2) 

lc 

where I (t) = return stroke current (Sunde’s model) 
ic 4 

a = 10 seconds-1 

b = 0.5 x 10 6 

10 = proportional to amount of charge 
displaced during return stroke 


Sunde’s model consists of two exponential terms and is relatively simple 
compared to some of the more recent models [pp 330-334 Uman, 1987]. For 
the time frames of interest in this study (t > 100 msec), it included 
the necessary temporal information required for late time transient 
electromagnetic analysis. 

The charge neutralized by the lightning current (charge 

perturbation) was described in terms of the temporal behavior of this 

current as follows (I. = dCL/dt): 

ic f 


Q f (t) 


t 


( 


I. (r)dx 
ic 

'0 


(3) 


where I. (t) = lightning current 

0 (t) = total amount of charge neutralized 

by lightning 
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The neutralization of charge by lightning occurs over the entire 
charge center and the possible effects of its spatial distribution must 
be considered in any thunderstorm electrical modeling. Its spatial 
profile will obviously strongly influence field behavior within the 
cloud, but at distances much larger than the charge centers extent 
(maximum feature size), its distribution will have little effect on the 
transient signatures [Baum, 1980]. Since the electric fields of 

interest were at least 20 km from the charge center (charge center’s 
maximum feature size - 2 km) there was a certain degree of freedom in 
the specification of the distribution of neutralized charge. Also, 
since there is virtually no published data describing the charge 
neutralization within the cloud the selection was even more arbitrary. 
A modified Gaussian spatial profile was chosen for the modeling and has 
been used in the modeling of many man-made and naturally occurring 
forced charge events [Baum, 1980]. 

The Gaussian profile used to describe the charge neutralized by 
lighting was also assumed for the charging caused by the thundercloud 
generator current. Here too, many possible profiles could have been 
used and for simplicity the Gaussian profile was selected. 

The mathematical description used in the modeling is shown below: 


±dp r (r, z,t)/<3t = ±(5Q (t)/3t ±oc)(f(r,z)) 

I i 

f(r.z) = (exp(-R Z /(2A)))/(2TrA) L5 
2 

where A = variance (A = 4000 m for simulations) 


(4) 


2 2 2 
R = r + (z - z’) 

z’ = altitude of charge center (6 or 10 km) 

& = term arising from steady state charging current 

(±1 A used in the modeling) 

A cylindrical coordinate system was used in the model with symmetry 
assumed about the vertical (z) axis (Figure 1). The effect of the 
magnitudes of the charge sources on the solutions will be discussed in 

a later section. 
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It is well known that the electric field recovery following 
lightning is a result of the neutralization of charge by the return 
stroke [pp 171-172 Uman, 1987 ]. The field's behavior is dictated 
primarily by the atmosphere's conductivity, the amount of charge 
neutralized by lightning and the magnitude of the thunderstorm’s 
current generator. 

Geometry of the Region 

Before describing the geometry, consider again the phenomenon cf 
interest: charge sources, located at altitudes no greater than 

approximately 10 km, induce electric fields throughout the atmosphere, 
but only those fields induced within the middle atmosphere were of 

interest in this study. Therefore, the model constructed must meet two 

criteria: 

1) If not obviously constrained, the geometrical limits of the 
model will approximate the entire atmosphere’s electrical effect on the 
regions where the simulations take place. 

2) The boundary conditions of the region will be electrically 

equivalent to those of the atmosphere. 

The region selected (Figure 1) is contained within a right circular 

cylinder with a radius of 80 km and height of 110 km. A discussion of 

how each of the boundaries was arrived at follows: 

- Lower Plate - The earth’s surface was electrically modeled as a 
perfect conductor. This assumption was based on the very large 

difference that exists between the earth’s conductivity and adjacent 

-3 -2 

atmosphere’s conductivity. Typically, values of 10 to 10 mhos/meter 

-14 -13 

are used for the earth’s conductivity, while 10 to 10 mhos/meter 

is the usual range of the adjacent atmosphere’s conductivity [pg 225, 

Voiiand, 1984]. This difference of more than 11 orders of magnitude 

makes the earth appear (electrically) as a perfect conductor. This 

assumption is commonly used in practical antenna engineering for 

-9 o 

frequencies whose range would correspond to time scales of 10 co 10 
seconds [Baum, 1980], ranges well beyond the maximum and minimum time 

scales of this research (minimum time step used in the research is 5 x 

-S 

10 seconds with a total duration 20 seconds.). 
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- Upper Plate - The selection of 110 km for the height of the 

upper boundary was a necessary consequence of the atmospheric 
conductivity structure being complicated by the Hall and Pederson 

components above an altitude of approximately 70 km [Volland, 1984]. 
The tensor conductivity components result when the mean free path and 

velocities of the charged carriers are sufficient to allow their 
trajectories to be altered by the effect of the earth’s magnetic field 
(Stratton, 1941], The finite element routine employed in the solution 
is not capable of modeling tensor conductivities, therefore an altitude 
limit must be set in the vicinity of 70 km. 

The UO-km altitude was selected based on the following 

considerations: 

1) The relative magnitude of the two tensor conductivity 

components is approximately proportional to the additional distance 
in altitude (beyond 70 km) considered. The maximum value of either of 

these components with respect to the parallel conductivity s magnitude 
(for the range of altitudes considered) is less than 20 percent 
[Volland, 1984]. 

2) The middle and upper atmospheric electric fields resulting from 

lightning (with the exception of the ~ 10' 4 propagating component) are 
approximately vertically oriented [Volland, 1984], i.e., the horizontal 

component is negligible. 

3) The off-diagonal tensor components of the conductivity will 
only interact with electric fields that are not aligned with the earth ^ 
magnetic field [Volland, 1984], Since the earth’s magnetic field with 
the exception of the equatorial regions, is primarily vertically aligned 
[Chalmers, 1967], the influence of both the Pederson and Hall components 
on the lightning induced vertical electric fields will be, to first 
order, negligible. 

An obvious concern is the influence this 110-km altitude limit may 
have on the simulations. To investigate the maximum possible error 

(assuming an electrically passive ionosphere) that this would introduce 

in the solutions, two sets of simulations were done with the 110-km 
upper plate electrically described by : 1) The vertical electric field 

and charge density are set to a value of zero. 2) The divergence of the 
electric field is set to the value of the charge density divided by the 
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and charge density are set to a value of zero. 2) The divergence of the 
electric field is set to the value of the charge density divided by the 
permittivity of free space (V*E = p/c). When the results were compared, 
little, if any, difference could be seen for the time frames of interest 
in any of the cases (all cases were subject to this test). 

The probable reason for this behavior is that, in general, for 
lightning-induced transients, the electrical properties of the 
atmosphere below the point of observation of the field, rather than 
above, govern the transient response [Hale and Baginski, 1987]. This 
may be explained by simply considering the fact the conductivity rapidly 
increases with altitude (i.e., resistivity decreasing), and therefore 
its influence (restrictive effect) on total global charge movement 
decreases. Hence, it seems reasonable to assume the middle atmosphere's 
simulated response to low-altitude charge perturbations is governed by 
the adjacent and lower altitude conductivity values. 

The conductivity used in the previously referenced high altitude 
simulations (Baginski et a l. 1988] was slightly modified to depict the 
observed conductivity [Blakeslee et al. y 1988] and used for the study in 
the form shown below (Figure 2): 

<r(z) = 1.3875 x 10~ 14 /(F1 + F2 + F3 + F4) 

where FI = 2.94/(exp(4.5 x 10 °x z)) 

F2 = 1.39/(exp(0.375 x 10~ 3 x z)) 

F3 = 0.369/(exp(0. 12987 x 10‘ 3 x z)) 

F4 = 9.0/(exp(2.777x 10' 3 x z)) 

Outer Cylindrical Surface - The lateral boundary had no distance 
constraint and could have been extended indefinitely. However, there 
exists a trade-off between accuracy and resolution: the smaller the 

model dimensions, the greater the accuracy in solving the differential 
equations. Therefore, the errors resulting from the adoption of finite 

boundaries for the model must be weighed against those resulting from 
degrading the numerical resolution of the simulation by involving too 
large a volume. 
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The simulations were found to be insensitive to increases in the 
radial limit beyond 70 km for all cases. No visible difference could be 
detected in the responses using either 70-km or 80-km radial boundaries 
when plotted together. Therefore, selecting a 80-km radial limit is a 
measure taken to provide additional confidence in the simulations. 


- Axis of Symmetry (z-axis) - Since r - 0 defines an axis of 
symmetry and since there are no discontinuities in the charge 

distribution, the derivative of the vertical electric field {dE^/dr) 
with respect to radial distance reduces to zero on the axis. 

The differential forms of the four boundaries are summarized as 
f ollows: 

1) At ARC = + 1 (z-axis), dE /dr - 0 

z 

2) At ARC = +2 (lower boundary), E^ = 0, equal potential surface 

3) At ARC = -3 (outer radial boundary), p = E^ = 0 

4) At ARC = +4 (upper boundary) V-E = p/e Q 

Maxwell’s Equations 

From Maxwell’s Equations a single equation is developed in which 
the electric field is dependent on the source charge and current 
densities as follows (Holzworth and Chiu, 1982): 


2 2 

Vx7xE = -u (crdE/dt + di /St + c d E/d t) (5) 

*0 so 


7 2 2 

Vp/e = v E - u {<rd E/dt + dJ /dt + c d E/d t) (6) 

K o o so 


where J = source current density associated with high intensity 
s 

lightning current, was neglected in the simulations 
p = charge density 


The continuity equation is derived by taking the divergence of the 
Maxwell current density: 

0 = 7 ■ 7xH = V * (crE dE/dt +J ) (7) 

o s 
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( 8 ) 


0 = pcr/c + 7<r-E + dp/d t + G 
o s 


where \7*J^ = G^ = source °f charge perturbation 

(deposition of charge by the lightning current and 
thunderstorm current generator) 


The above equations were derived from Maxwell* s equations and 
were used to describe the electrodynamic response of the atmosphere to 
the thunderstorm.. It is important to note the primary difference in 
the forcing function (charge related) used here relative to forcing 

functions used for most spheric analysis (typically current-related 
forcing functions). A current related transient (e.g., current loop) is 
allowable in a charge neutral region, but charge movement through space 
does not necessarily imply a zero volume charge density. Furthermore, a 
charge perturbation may be considered separately from the currents which 
cause it (charge perturbation implies charge generation at a point in 

space). In the latter case charge motion is required (conservation of 

charge), but the relative contribution of each (charge versus current) 
mechanism to the resulting total electric field may be analyzed 

separately. Initially (prior to the first iterative time step of the 

transient analysis), no charge is displaced and the transient electric 
field everywhere is assumed zero. 

Equations 6 and 8 used in the modeling are linear partial 
differential equations. This fact greatly reduced the degree of 

complexity in the thunderstorm electrical modeling. Linearity allows 

the use of superposition to construct the total solution desired based 
on each source's individual response. The amplitude of the solution is 

also directly proportional to the source’s magnitude. This was 

especially important in the comparative analysis when the specific 
thunderstorm’s current generator and amount of charge neutralized by 
lightning were derived. 
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SIMULATIONS 


The vertical component of the simulated Maxwell current density 
signatures are shown in Figure 3 for altitudes of 46, 53, and 60 km at a 
radial distance of 15 km from the vertical axis of symmetry. The 
waveforms show only slight differences in both magnitude and temporal 

character. Vertical electric field signatures at the same positions 
show relatively large variations in both magnitude and 
temporal character (Figure 4). 

Figure 5 and 6 identify the normalized Maxwell current density and 
electric field for altitudes of 40-60 km and radial distances of 15-30 
km at t - 200 milliseconds. The interesting feature is that the Maxwell 

current density and electric field are not aligned. This observation 
is of importance to experimentalists who typically infer field 

characteristics from a single measured quantity (E and J ). 

The simulated vertical electric field signatures (Figure 4) are 
based on conductivity measurements obtained concurrently with with 
transient electric field data. There are two traits common to both 
measurements and simulations that are not explainable in terms of simple 
relaxation time solutions: the significant time delay prior to the onset 
of the maximum field strength, and the relatively long duration of the 
transient. This fact indicates that the use of relaxation time 
solutions would not be valid for the study conducted. 

CONCLUSIONS 

The major contribution of the study Is to emphasis the general 
importance modeling has in identifying the allowable atmospheric 

transient ^ behavior. This reduces the researchers need to rely on 
assumptions that are in many cases not rigorously justifiable. It is 

probable that the next major advances in atmospheric electricity will 
involve numerical simulation predictions of otherwise unexpected 

behavior. 

Lightning researchers may benefit from the relative spatial and 
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temporal invariance of the Maxwell current density with respect to the 
corresponding electric field signature. It is likely that this type of 
characterization would provide an alternative fundamental description of 
the lightning event that would greatly reduce the complexity of much 
current thunderstorm electrical research. 
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Figure 2 Composite Conductivity Profile Based on Measured Data 
Obtained August 8, 1981 at Wallops Island, Virginia. 
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Figure 3 Maxwell Current Density Simulations at a 
Radial Distance of 15 km for Altitudes of 
46, 53, and 60 km. 
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Figure 4 Electric Field Simulations at a Radial Distance 
of 15 km for Altitudes of 46, 53, and 60 km. 
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Figure 5 Normalized Simulated Electric Field Mapping for 
t * 200 milliseconds . 
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